library(forecast)
library(fpp2)
## Loading required package: fma
##
## Attaching package: 'fma'
## The following object is masked from 'package:plyr':
##
## ozone
## Loading required package: expsmooth
library(tidyverse)
library(tseries)
detach("package:dplyr", unload=TRUE)
## Warning: 'dplyr' namespace cannot be unloaded:
## namespace 'dplyr' is imported by 'tigris', 'broom', 'tidycensus' so cannot be unloaded
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:xts':
##
## first, last
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(scales)
##
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
##
## discard
## The following object is masked from 'package:readr':
##
## col_factor
library(GGally)
##
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
##
## nasa
## The following object is masked from 'package:fma':
##
## pigs
library(e1071)
library(fpp2)
theme_set(theme_minimal())
3.1 Some simple forecasting methods
Average method
\[\hat{y}_{T+h|T} = (y_{1} + \dots + y_{T})/T\]
Naive method
\[\hat{y}_{T+h|T} = y_{t}\]
Seasonal Naive method
\[\hat{y}_{T+h|T} = y_{T+h-m(k+1)}\]
For example, with monthly data, the forecast for all future February values is equal to the last observed February value.
Drifted method
\[\hat{y}_{T+h|T} = y_T + \frac{h}{T-1}\sum_{t=2}^T(y_t-y_{t-1}) = y_T +h \frac{y_T - y_1}{T-1}\] This is equivalent to drawing a line between the first and last observations, and extrapolating it into the future.
x <- ts(1:20, start = 1990, frequency = 4)
autoplot(x)
autoplot(naive(x, 5))
autoplot(rwf(x, 4))
autoplot(snaive(x, 4))
beer2 <- window(ausbeer, start = 1992, end = c(2007, 4))
autoplot(beer2) +
autolayer(meanf(beer2, h = 11),
series = "Mean", PI = F)+
autolayer(naive(beer2, h = 11),
series = "Naive", PI = F) +
autolayer(snaive(beer2, h = 11),
series = "Seasonal Naive", PI = F) +
autolayer(rwf(beer2, h = 11, drift = T),
series = "Naive with Drift", PI = F) +
labs(x = "Year", y = "Megalitres",
title = "Foecasts for quarterly beer production")
These simple method will be consider as benchmarks. all new methods will be compare to these mmethod, to see whether method are better or not.
3.2 Transformation and adjustments
Calender adjustments
autoplot(milk)
dframe <- cbind(monthly = milk,
DailyAverage = milk / monthdays(milk))
autoplot(dframe, facet = T) +
labs(x = "year", y = "Pounds",
title = "Milk production per cow")
Population adjustments
Inflation Adjustments
Mathematical Transformations
log transformations, power transformations, Box Cox transformations
lambda <- BoxCox.lambda(elec)
lambda
## [1] 0.2654076
autoplot(BoxCox(elec,lambda))
kpss.test(elec)
## Warning in kpss.test(elec): p-value smaller than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: elec
## KPSS Level = 7.9424, Truncation lag parameter = 5, p-value = 0.01
dframe <- cbind(month = elec,
DailyAverage = elec / monthdays(elec))
autoplot(dframe, facets = T) +
labs(x = "year")
Bias adjustments
fc <- rwf(eggs, drift = T, lambda = 0, h = 50, level = 80)
fc2 <- rwf(eggs, drift = T, lambda = 0, h = 50, level = 80, biasadj = T)
autoplot(eggs) +
autolayer(fc, series = "Simple back transformation") +
autolayer(fc2, series = "Bias adjusted", PI = F) +
guides(color = guide_legend(title = "Forecast"))
Check residual normality like linear regression does 1. residual uncorrelated 2. Residuals has a zero mean 3. Residuals have constant variance 4. The residuals are normally distributed.
Example: Forecasting the Google dialy Closing stock price.
\[e_{t} = y_t - \hat{y_t} = y_t - y_{t-1}\]
autoplot(goog200) +
labs(x = "Day", y = "Close Price (US$)",
title = "Google Stock (daily ending 6 December 2013)")
we can chekc the residual obtain from forecasting this series using the naive method
# plot the naive forecasting
goog_naive <- naive(goog200, 100)
autoplot(goog200) +
autolayer(goog_naive, se = F, series = "Naive") +
labs(x = "Day", y = "Close Price (US$)",
title = "Google Stock (daily ending 6 December 2013)")
plot the residuals
res <- residuals(goog_naive)
t.test(res)
##
## One Sample t-test
##
## data: res
## t = 1.5892, df = 198, p-value = 0.1136
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.1678207 1.5612704
## sample estimates:
## mean of x
## 0.6967249
autoplot(res) +
labs(x = "Day", y = "",
title = "Residuals from naive method")
gghistogram(res) +
labs(title = "Histogram of residuals") +
stat_function(fun = dnorm,
color = "red",
args = list(mean = mean(res, na.rm = T),
sd = sd(res, na.rm = T)))
## Warning: Removed 1 rows containing non-finite values (stat_bin).
ggAcf(res) +
labs(title = "ACF of residuals")
These graphs show that the naive method produces forecasts that appear to account for all avaiable info. The mean of the residuals is calose to zero, in t test result. There is no correlation in the residual series, in acf test no value exceed \[\pm2/\sqrt{n}\], n is the length of the Time Series. The times plot of the residuals show that the variation of teh residuals stay much the same across the historical data, aprat from the outlier, and therefore residual can be treated as constant.
Consequently, forecasts from naive method is a good fit, but the prediction interval that are computed assuming a normal distributuon may be inaccurate.
# lag = h and fitdf = K
Box.test(res, lag = 10, fitdf = 0)
##
## Box-Pierce test
##
## data: res
## X-squared = 10.611, df = 10, p-value = 0.3886
Box.test(res, lag = 10, fitdf = 0, type = "Lj")
##
## Box-Ljung test
##
## data: res
## X-squared = 11.031, df = 10, p-value = 0.3551
Since both Q and \(Q^*\), the result are no significatnt. Thus, we can clude that the residuals are not distriguishable from a white noise sieres.
All those methods for checking residuals are conveniently in R function checkresiduals()
checkresiduals(naive(goog200))
##
## Ljung-Box test
##
## data: Residuals from Naive method
## Q* = 11.031, df = 10, p-value = 0.3551
##
## Model df: 0. Total lags used: 10
start(ausbeer)
## [1] 1956 1
end(ausbeer)
## [1] 2010 2
window(ausbeer, start = 1995)
## Qtr1 Qtr2 Qtr3 Qtr4
## 1995 426 408 416 520
## 1996 409 398 398 507
## 1997 432 398 406 526
## 1998 428 397 403 517
## 1999 435 383 424 521
## 2000 421 402 414 500
## 2001 451 380 416 492
## 2002 428 408 406 506
## 2003 435 380 421 490
## 2004 435 390 412 454
## 2005 416 403 408 482
## 2006 438 386 405 491
## 2007 427 383 394 473
## 2008 420 390 410 488
## 2009 415 398 419 488
## 2010 414 374
window(ausbeer, start = length(ausbeer) - 4*5)
## Warning in window.default(x, ...): 'start' value not changed
## Qtr1 Qtr2 Qtr3 Qtr4
## 1956 284 213 227 308
## 1957 262 228 236 320
## 1958 272 233 237 313
## 1959 261 227 250 314
## 1960 286 227 260 311
## 1961 295 233 257 339
## 1962 279 250 270 346
## 1963 294 255 278 363
## 1964 313 273 300 370
## 1965 331 288 306 386
## 1966 335 288 308 402
## 1967 353 316 325 405
## 1968 393 319 327 442
## 1969 383 332 361 446
## 1970 387 357 374 466
## 1971 410 370 379 487
## 1972 419 378 393 506
## 1973 458 387 427 565
## 1974 465 445 450 556
## 1975 500 452 435 554
## 1976 510 433 453 548
## 1977 486 453 457 566
## 1978 515 464 431 588
## 1979 503 443 448 555
## 1980 513 427 473 526
## 1981 548 440 469 575
## 1982 493 433 480 576
## 1983 475 405 435 535
## 1984 453 430 417 552
## 1985 464 417 423 554
## 1986 459 428 429 534
## 1987 481 416 440 538
## 1988 474 440 447 598
## 1989 467 439 446 567
## 1990 485 441 429 599
## 1991 464 424 436 574
## 1992 443 410 420 532
## 1993 433 421 410 512
## 1994 449 381 423 531
## 1995 426 408 416 520
## 1996 409 398 398 507
## 1997 432 398 406 526
## 1998 428 397 403 517
## 1999 435 383 424 521
## 2000 421 402 414 500
## 2001 451 380 416 492
## 2002 428 408 406 506
## 2003 435 380 421 490
## 2004 435 390 412 454
## 2005 416 403 408 482
## 2006 438 386 405 491
## 2007 427 383 394 473
## 2008 420 390 410 488
## 2009 415 398 419 488
## 2010 414 374
tail(ausbeer, 4*5)
## Qtr1 Qtr2 Qtr3 Qtr4
## 2005 408 482
## 2006 438 386 405 491
## 2007 427 383 394 473
## 2008 420 390 410 488
## 2009 415 398 419 488
## 2010 414 374
Forecast error is the difference between an observation and its foecast
\[e_{T+h} = y_{T+h} - \hat{y}_{T+h|T}\] forecast error different from residuals in two ways
we can measure forecast accuracy by summarising the forecast erro in different ways.
Mean absolute error: MAE = mean\(|e_t|\) Root mean squared error: RMSE = \(\sqrt{mean(e_t^2)}\)
When comparing forecast methods applied to a single time series, or to serveral time series with the same units, the MAE is popular as it it easy to both understand and compute. A forecast method that minise the MAE will lead to forecasts of the median, while minimising the RMSE will lead to forecast of the mean. Consequently, the RMSE is also widely used, depsite being more difficult to interpret.
\(p_t = 100e_t/y_t\), mean unit-free
Mean absolute percentage error: MAPE = mean(\(|p_t|\)) Potential drawdown, being infinite or undefined if \(y_t = 0\) for any t in the period of interest.
Scaled errors as an alternative to using percentage errors when comparing forecast accuracy across series with different units.
Seaonal dataset
beer2 <- window(ausbeer,start=1992,end=c(2007,4))
beerfit1 <- meanf(beer2,h=10)
beerfit2 <- rwf(beer2,h=10)
beerfit3 <- snaive(beer2,h=10)
autoplot(window(ausbeer, start=1992)) +
autolayer(beerfit1, series="Mean", PI=FALSE) +
autolayer(beerfit2, series="Naïve", PI=FALSE) +
autolayer(beerfit3, series="Seasonal naïve", PI=FALSE) +
labs(x = "year", y = "Megalitres", title = "Forecasts for quarterly beer production")
beer3 <- window(ausbeer, start = 2008)
accuracy(beerfit1, beer3)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.000 43.62858 35.23438 -0.9365102 7.886776 2.463942
## Test set -13.775 38.44724 34.82500 -3.9698659 8.283390 2.435315
## ACF1 Theil's U
## Training set -0.10915105 NA
## Test set -0.06905715 0.801254
accuracy(beerfit2, beer3)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.4761905 65.31511 54.73016 -0.9162496 12.16415 3.827284
## Test set -51.4000000 62.69290 57.40000 -12.9549160 14.18442 4.013986
## ACF1 Theil's U
## Training set -0.24098292 NA
## Test set -0.06905715 1.254009
accuracy(beerfit3, beer3)
## ME RMSE MAE MPE MAPE MASE
## Training set -2.133333 16.78193 14.3 -0.5537713 3.313685 1.0000000
## Test set 5.200000 14.31084 13.4 1.1475536 3.168503 0.9370629
## ACF1 Theil's U
## Training set -0.2876333 NA
## Test set 0.1318407 0.298728
Base on the accuracy result, we can see the thir mothod, Seaonal naive method is the most accuracy.
non-seasonal dataset
# three different method
googfc1 <- meanf(goog200, h = 40)
googfc2 <- rwf(goog200, h = 40)
googfc3 <- rwf(goog200, h = 40, drift = T)
autoplot(subset(goog, end = 240)) +
autolayer(googfc1, PI = F, series = "Mean") +
autolayer(googfc2, PI = F, series = "Naive") +
autolayer(googfc3, PI = F, series = "Drift") +
labs(x = "day", y = "Close price (US$)",
title = "Google Stock price forecast")
Find accuracy
googtest <- window(goog, start = 201, end = 240)
accuracy(googfc1, googtest)
## ME RMSE MAE MPE MAPE
## Training set -4.296286e-15 36.91961 26.86941 -0.6596884 5.95376
## Test set 1.132697e+02 114.21375 113.26971 20.3222979 20.32230
## MASE ACF1 Theil's U
## Training set 7.182995 0.9668981 NA
## Test set 30.280376 0.8104340 13.92142
accuracy(googfc2, googtest)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.6967249 6.208148 3.740697 0.1426616 0.8437137 1.000000
## Test set 24.3677328 28.434837 24.593517 4.3171356 4.3599811 6.574582
## ACF1 Theil's U
## Training set -0.06038617 NA
## Test set 0.81043397 3.451903
accuracy(googfc3, googtest)
## ME RMSE MAE MPE MAPE
## Training set -5.998536e-15 6.168928 3.824406 -0.01570676 0.8630093
## Test set 1.008487e+01 14.077291 11.667241 1.77566103 2.0700918
## MASE ACF1 Theil's U
## Training set 1.022378 -0.06038617 NA
## Test set 3.119002 0.64732736 1.709275
The best method, which have the lowest value in all summarize error term method, is drift method.
e <- tsCV(goog200, rwf, drift = T, h = 1)
# RMSE: Root mean squared error
sqrt(mean(e^2, na.rm = T))
## [1] 6.233245
sqrt(mean(residuals(rwf(goog200, drift = 1))^2, na.rm = T))
## [1] 6.168928
The RMSE from the residuals is smaller, as the corrsponding “forecast” are based on a model fitted to the entire data set, rather than being true forecasts.
A good way to choose the best forecasting model is to find the model with smallest RMSE computed using time series cross-validation.
error <- goog200 %>%
tsCV(forecastfunction = rwf, drift = T, h = 1)
error^2 %>%
mean(na.rm = T) %>%
sqrt()
## [1] 6.233245
error <- goog200 %>%
rwf(drift = T) %>%
residuals()
error^2 %>%
mean(na.rm = T) %>%
sqrt()
## [1] 6.168928
error <- tsCV(goog200, forecastfunction = naive, h = 8)
# compute the MSE value and remove missing values
mse <- colMeans(error^2, na.rm = T)
data.frame(h = 1:8, MSE = mse) %>%
ggplot(aes(x = h, y = MSE)) +
geom_point()
# rmse
e <- tsCV(goog200, rwf, drift = T, h = 8)
rmse <- sqrt(colMeans(error^2, na.rm = T))
rmse
## h=1 h=2 h=3 h=4 h=5 h=6 h=7
## 6.208148 8.578760 10.730162 12.845496 14.655090 16.082789 17.510978
## h=8
## 19.150658
data.frame(h = 1:8, RESE = rmse) %>%
ggplot(aes(x = h, y = RESE)) +
geom_point()
\(\hat{y}_{T+h|T}\pm{c}\hat{\sigma}_h\)
# naive method residuals
residuals <- residuals(naive(goog200))
sd(residuals, na.rm = T)
## [1] 6.184487
Common feature of predication intervals is that they increase in length as the forecast horizon increases. The future ahead we forecast, the more uncertainity is associated with the forecast, and thus the wider the predication intervals.
Four benchmark methods
Note that when h = 1 and T is large, means you need to only predict one value and you have large dataset like 200, these all give the same approximate value \(\hat{\sigma}\)
naive(goog200, h = 10)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 201 531.4783 523.5222 539.4343 519.3105 543.6460
## 202 531.4783 520.2267 542.7298 514.2705 548.6861
## 203 531.4783 517.6980 545.2586 510.4031 552.5534
## 204 531.4783 515.5661 547.3904 507.1428 555.8138
## 205 531.4783 513.6880 549.2686 504.2704 558.6862
## 206 531.4783 511.9900 550.9666 501.6735 561.2830
## 207 531.4783 510.4285 552.5280 499.2854 563.6711
## 208 531.4783 508.9751 553.9814 497.0627 565.8939
## 209 531.4783 507.6101 555.3465 494.9750 567.9815
## 210 531.4783 506.3190 556.6375 493.0005 569.9561
goog200[200] + sd(residuals(naive(goog200)), na.rm = T) * qnorm(0.1, lower.tail = F)
## [1] 539.404
autoplot(naive(goog200))
### Predication intervals from bootstrpped residuals
When a normal distribution for the forecast errors is an unreasonable assumption, one alternative is to use bootstrapping, which only assumes that the forecast errors are uncorrelated
Forecast error: \(e_t = y_t - \hat{y}_{t|t-1}\) which can re-rwite this as \(y_t = \hat{y}_{t|t-1} + e_t\)
==================Review
naive(goog200, bootstrap = T)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 201 531.4783 526.0076 537.0061 523.2307 541.0597
## 202 531.4783 523.1504 539.8016 519.8481 546.2497
## 203 531.4783 521.0374 541.6380 517.0867 550.3955
## 204 531.4783 519.1055 543.4188 514.4375 554.0778
## 205 531.4783 517.6330 544.9533 511.8730 563.3855
## 206 531.4783 516.3541 546.7031 509.7720 581.3736
## 207 531.4783 514.8630 547.9053 507.8135 581.9491
## 208 531.4783 513.7869 549.3288 505.9135 586.0329
## 209 531.4783 512.4148 550.4386 503.6003 586.3873
## 210 531.4783 510.9900 551.6743 502.1216 588.5228
In this case, they are similary (but not identical) to the predication intervals based on teh normal distribution.
If a transfromation has been used, then the predication interval should be computed on the transformed scale, and the end points back-transformation to give a predication interval on teh original scale.
forecast(ausbeer, h = 4)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2010 Q3 404.6025 385.8695 423.3355 375.9529 433.2521
## 2010 Q4 480.3982 457.5283 503.2682 445.4216 515.3748
## 2011 Q1 417.0367 396.5112 437.5622 385.6456 448.4277
## 2011 Q2 383.0996 363.5063 402.6930 353.1341 413.0651
# US net electricity generation
autoplot(usnetelec)
lambda <- BoxCox.lambda(usnetelec)
autoplot(BoxCox(usnetelec, lambda))
kpss.test(BoxCox(usnetelec, lambda))
## Warning in kpss.test(BoxCox(usnetelec, lambda)): p-value smaller than
## printed p-value
##
## KPSS Test for Level Stationarity
##
## data: BoxCox(usnetelec, lambda)
## KPSS Level = 1.4583, Truncation lag parameter = 3, p-value = 0.01
# Quarterly US GDP
autoplot(usgdp)
lambda <- BoxCox.lambda(usgdp)
autoplot(BoxCox(usgdp, lambda))
kpss.test(BoxCox(usgdp, lambda))
## Warning in kpss.test(BoxCox(usgdp, lambda)): p-value smaller than printed
## p-value
##
## KPSS Test for Level Stationarity
##
## data: BoxCox(usgdp, lambda)
## KPSS Level = 4.8114, Truncation lag parameter = 4, p-value = 0.01
# monthly copper prices
autoplot(mcopper)
lambda <- BoxCox.lambda(mcopper)
autoplot(BoxCox(mcopper, lambda))
kpss.test(BoxCox(mcopper, lambda))
## Warning in kpss.test(BoxCox(mcopper, lambda)): p-value smaller than printed
## p-value
##
## KPSS Test for Level Stationarity
##
## data: BoxCox(mcopper, lambda)
## KPSS Level = 6.2659, Truncation lag parameter = 6, p-value = 0.01
# monthly US domestic enplanements
autoplot(enplanements)
lambda <- BoxCox.lambda(enplanements)
autoplot(BoxCox(enplanements, lambda))
autoplot(cangas)
lambda <- BoxCox.lambda(cangas)
autoplot(BoxCox(cangas, lambda))
kpss.test(BoxCox(cangas, lambda))
## Warning in kpss.test(BoxCox(cangas, lambda)): p-value smaller than printed
## p-value
##
## KPSS Test for Level Stationarity
##
## data: BoxCox(cangas, lambda)
## KPSS Level = 7.2459, Truncation lag parameter = 6, p-value = 0.01
Because the variance is not simple increase or decrease with series, it jump during 70-90 and decrease afterward.
retail_data <- readxl::read_excel("/Users/yifeiliu/Documents/R/data/book_exercise/forecasting/retail.xlsx", skip = 1)
myts <- ts(retail_data[, "A3349873A"],
frequency = 12, start = c(1982, 4))
autoplot(myts)
lambda <- BoxCox.lambda(myts)
autoplot(BoxCox(myts, lambda))
kpss.test(BoxCox(myts, lambda))
## Warning in kpss.test(BoxCox(myts, lambda)): p-value smaller than printed p-
## value
##
## KPSS Test for Level Stationarity
##
## data: BoxCox(myts, lambda)
## KPSS Level = 5.8234, Truncation lag parameter = 5, p-value = 0.01
# unemployment benefits in Australia
autoplot(dole)
dole_cal <- dole / monthdays(dole)
autoplot(dole_cal)
# transformation adjust for calender, population, inflation seems appropriate
# monthly accident deaths in USA
autoplot(usdeaths)
usdeaths_cal <- cbind(usdeath = usdeaths,
calender = usdeaths / monthdays(usdeaths))
autoplot(usdeaths_cal, facets = T)
# calender transformation seem reasonable, if we have further dataset such as ppulation. maybe population tranform is also reasonable
# quarterly clay brick production
autoplot(bricksq)
lambda <- BoxCox.lambda(bricksq)
autoplot(BoxCox(bricksq, lambda))
# all the transfromation and adjustment I know don't seem appropriate to apply on this dataset.
beer <- window(ausbeer, start = 1992)
fc <- snaive(beer)
autoplot(fc)
res <- residuals(fc)
autoplot(res)
checkresiduals(res)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.
t.test(res)
##
## One Sample t-test
##
## data: res
## t = -0.81286, df = 69, p-value = 0.4191
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -5.428075 2.285218
## sample estimates:
## mean of x
## -1.571429
To test if the residuals are white noise we can do Portmanteau test for autocorrelation
# lag = h, h = 2*m, m = 4. fitdf = k, k = parameters in the model, model is snave, so k = 1
Box.test(res, lag = 2*4, fitdf = 1, type = "Lj")
##
## Box-Ljung test
##
## data: res
## X-squared = 32.269, df = 7, p-value = 3.621e-05
The p value is relative large, which mean the result are signiciant, so we cannot conclude that the residuals are distringuishable from a white noise series.
# WWWusage internet usage per minute
autoplot(WWWusage)
usage_mean <- rwf(WWWusage)
mean_res <- residuals(usage_mean)
checkresiduals(mean_res)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.
usage_snaive <- snaive(WWWusage)
snaive_res <- residuals(usage_mean)
checkresiduals(snaive_res)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.
usage_naive <- naive(WWWusage)
naive_res <- residuals(usage_naive)
checkresiduals(usage_naive)
##
## Ljung-Box test
##
## data: Residuals from Naive method
## Q* = 145.58, df = 10, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 10
# both naive or snaive is not approporate
# bricksq
autoplot(bricksq)
bricksq_naive <- naive(bricksq)
naive_res <- residuals(bricksq_naive)
checkresiduals(naive_res)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.
bricksq_snaive <- snaive(bricksq)
snaive_res <- residuals(bricksq_snaive)
checkresiduals(snaive_res)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.
# both are not appropriated
myts.train <- window(myts, end=c(2010,12))
myts.test <- window(myts, start=2011)
autoplot(myts) +
autolayer(myts.train, series="Training") +
autolayer(myts.test, series="Test")
fc <- snaive(myts.train)
accuracy(fc,myts.test)
## ME RMSE MAE MPE MAPE MASE
## Training set 7.772973 20.24576 15.95676 4.702754 8.109777 1.000000
## Test set 55.300000 71.44309 55.78333 14.900996 15.082019 3.495907
## ACF1 Theil's U
## Training set 0.7385090 NA
## Test set 0.5315239 1.297866
checkresiduals(fc)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 624.45, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24
kurtosis(residuals(fc), na.rm = T)
## [1] 1.311521
# we can see the residuals apprea to be autocorrelated and distribution also have excess kurtosis.
# use cross validation to test accuracy
e <- tsCV(myts, snaive, h = 4)
sqrt(mean(e^2, na.rm = T))
## [1] 25.50797
resids <- myts %>%
snaive() %>%
residuals()
resids^2 %>%
mean(na.rm = T) %>%
sqrt()
## [1] 25.68062
# test sensitive to acuracy measurement
percent(length(myts.test) / length(myts))
## [1] "9.45%"
test sensitive to acuracy measurement, we can see the tes data set we only use 9.45% of original dataset. Convention of train/test should be .8/.2. So we can these how snaive model perform in these split case
myts.train <- window(myts, end = c(2007, 6))
myts.test <- window(myts, start = c(2007, 7))
autoplot(myts) +
autolayer(myts.train, series="Training") +
autolayer(myts.test, series="Test")
fc <- snaive(myts.train)
accuracy(fc,myts.test)
## ME RMSE MAE MPE MAPE MASE
## Training set 8.969072 19.52551 15.28729 5.450456 8.272556 1.000000
## Test set 19.358333 24.14890 20.78333 6.124900 6.664263 1.359518
## ACF1 Theil's U
## Training set 0.7198404 NA
## Test set 0.6031662 0.4537988
checkresiduals(fc)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 584.22, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24
kurtosis(residuals(fc), na.rm = T)
## [1] 1.982122
20/80 split produce accuracy result better than 10/90 split
# quarterly visitor nights for various region of australia
autoplot(visnights)
qldmetro <- (visnights[, "QLDMetro"])
qldmetro_tr1 <- window(qldmetro, end = c(2015, 4))
qldmetro_tr2 <- window(qldmetro, end = c(2014, 4))
qldmetro_tr3 <- window(qldmetro, end = c(2013, 4))
fc1 <- snaive(qldmetro_tr1)
fc2 <- snaive(qldmetro_tr2)
fc3 <- snaive(qldmetro_tr3)
accuracy(fc1, qldmetro)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.02006107 1.0462821 0.8475553 -0.2237701 7.976760 1.0000000
## Test set 0.56983879 0.9358727 0.7094002 4.6191866 6.159821 0.8369957
## ACF1 Theil's U
## Training set 0.06014484 NA
## Test set 0.09003153 0.4842061
accuracy(fc2, qldmetro)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0161876 1.0735582 0.8809432 -0.2744747 8.284216 1.0000000
## Test set 0.3669560 0.8516202 0.6644561 2.8431247 6.085501 0.7542553
## ACF1 Theil's U
## Training set 0.066108793 NA
## Test set -0.002021023 0.5102527
accuracy(fc3, qldmetro)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.007455407 1.074544 0.8821694 -0.5625865 8.271365 1.0000000
## Test set 0.411850940 1.035878 0.8800104 4.4185560 8.635235 0.9975526
## ACF1 Theil's U
## Training set 0.07317746 NA
## Test set -0.21876890 0.8365723
we can see the result from training dataset 2 produce the lowest MAPE score.
autoplot(dowjones)
autoplot(rwf(dowjones, drift = T))
# this forecast are identical to extending the line drawn between the first and last observation
autoplot(rwf(dowjones, drift = T))